﻿using PGMHelper;
using System.Windows.Forms.DataVisualization.Charting;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Spectrum
{
    /// <summary>
    /// 人工波对象
    /// </summary>
    public class ArtiWaveObj
    {
        #region Variables

        /// <summary>
        /// 人工波名称
        /// </summary>
        protected string Name { set; get; }

        /// <summary>
        /// 平稳高斯过程
        /// </summary>
        private List<StationaryGaussianParam> StationaryGaussianList { set; get; }

        /// <summary>
        /// 控制周期
        /// </summary>
        private List<double> ControlPeriods { set; get; }

        /// <summary>
        /// 加速度时程
        /// </summary>
        public AccelSet accels { set; get; }

        /// <summary>
        /// 包络函数对象
        /// </summary>
        private EnvelopeFunction Ft { set; get; }

        /// <summary>
        /// 国家反应谱
        /// </summary>
        private SpectSet codeSpectrum { set; get; }

        /// <summary>
        /// 地震动反应谱
        /// </summary>
        public SpectSet waveSpectrum { set; get; }

        /// <summary>
        /// 事件间隔
        /// </summary>
        private double TimeInterval { set; get; }

        /// <summary>
        /// 地震动持时
        /// </summary>
        private int Duration { set; get; }

        /// <summary>
        /// 迭代次数
        /// </summary>
        private int IterationCount { set; get; }

        /// <summary>
        /// 时间编号最大值
        /// </summary>
        private int TimeIndexMax { set; get; }

        #endregion

        /// <summary>
        /// 设定非平稳随机过程
        /// </summary>
        private void SetAccels()
        {
            //获得时程
            this.accels = this.GetAccels();
            //获得峰值加速度(g)
            var pga = this.accels.OriginAccelPGA;
            //目标峰值加速度
            var targetPGA = 1 / this.codeSpectrum.GetReductionFactor();
            //对时程进行放缩
            this.StationaryGaussianList.ForEach(gauss => gauss.ModifyPowerSpectrum((float)(targetPGA / pga)));
            //获得时程
            this.accels = this.GetAccels();
            //获得峰值加速度(g)
            pga = this.accels.OriginAccelPGA;
        }

        /// <summary>
        /// 获得加速度时程
        /// </summary>
        /// <returns></returns>
        private AccelSet GetAccels()
        {
            //初始化激励
            var accelList = new List<AccelParam>((int)(this.Duration / TimeInterval + 1));
            accelList.Add(new AccelParam { Accel = 0, Time = 0 });
            //遍历
            for (int index = 1; index <= this.TimeIndexMax; index++)
            {
                //平稳高斯过程
                var Ast = new double();
                //获得时间
                var time = index * this.TimeInterval;
                //遍历列表累加
                this.StationaryGaussianList.ForEach(gauss => Ast += gauss.GetAs(index - 1));
                //转为非平稳随机过程
                accelList.Add(new AccelParam
                {
                    Accel = (float)(Ast * this.Ft.GetFt(time)),
                    Time = (float)time
                });
            }
            //设定加速度
            return new AccelSet(accelList);
        }

        /// <summary>
        /// 设定反应谱
        /// </summary>
        /// <returns></returns>
        private SpectSet SetSpectrum(bool isFinal)
        {
            //构造反应谱求解器
            var spectrumAlog = new SpectrumMethodParams(this.accels, (float)ArtiWaveParam.DampRatio,
                ArtiWaveParam.PeriodLength + 1, (float)ArtiWaveParam.PeriodDelta, (float)this.TimeInterval);
            //清空参数
            spectrumAlog.ClearNewMark();
            //求解反应谱
            this.waveSpectrum = isFinal? spectrumAlog.GetSpectParam(true, ArtiWaveParam.SpectrumDtMax):
                spectrumAlog.GetSpectParamFast(true, this.ControlPeriods);
            //反应谱调幅
            if (isFinal) this.waveSpectrum.SpectrumAmplitude((float)ArtiWaveParam.SpecturmModifiedFactor);
            //返回对象
            return this.waveSpectrum;
        }

        /// <summary>
        /// 设定高斯平稳过程
        /// </summary>
        private void SetStationaryGaussianBySameW()
        {
            //初始化
            this.StationaryGaussianList = new List<StationaryGaussianParam>();
            this.ControlPeriods = new List<double>();
            //构造随机数
            Random rd = new Random();
            //确定功率增量
            var deltaW = //2 * Math.PI / this.Duration;
               Math.Min(2 * Math.PI / this.Duration,
               2 * Math.PI / (ArtiWaveParam.PeriodLength) - 2 * Math.PI / (ArtiWaveParam.PeriodLength + ArtiWaveParam.DeltaPeriodCreateGaussian));
            //最大功率
            var Wmax = 2 * Math.PI / (ArtiWaveParam.MinDwFactor * ArtiWaveParam.ModifiedPeriodDelta);
            //遍历构造高斯平稳过程
            for (var w = 2 * Math.PI / (ArtiWaveParam.PeriodLength + 1); w <= Wmax; w += deltaW)
            {
                //获得基础相位角
                var testThetaK = 2 * Math.PI / 100 * rd.Next(100);
                //初始化三角级数
                var Gaussian = new StationaryGaussianParam(w, deltaW, testThetaK, this.Duration, this.codeSpectrum);
                //设定三角级数
                Gaussian.SetCosValues(this.TimeIndexMax, this.TimeInterval);
                this.StationaryGaussianList.Add(Gaussian);
                //计算周期值
                var period = 2 * Math.PI / w;
                //计算前后周期值
                var periodMin = (int)(period / ArtiWaveParam.ModifiedPeriodDelta) * ArtiWaveParam.ModifiedPeriodDelta;
                //确定目标周期值
                var preiodTarget = period - periodMin < 0.5 * ArtiWaveParam.ModifiedPeriodDelta ? periodMin : periodMin + ArtiWaveParam.ModifiedPeriodDelta;
                preiodTarget = Math.Round(preiodTarget, 2);
                //判断列表是否存在
                if (!this.ControlPeriods.Contains(preiodTarget)) this.ControlPeriods.Add(preiodTarget);
            }
            //逆序排列
            this.ControlPeriods.Reverse();
            this.StationaryGaussianList.Reverse();
        }

        /// <summary>
        /// 仅调整某目标频率下的分量
        /// </summary>
        private void ModifyPowerSpectrumA()
        {
            //周期增量
            var deltaPeiod = (float)ArtiWaveParam.ModifiedPeriodDelta;
            //高斯列表序号
            int GuassIndex = 0;
            //新的周期编号
            var nextGuassPeriod = 0.0;
            //遍历地震动反应谱
            foreach (float period in this.ControlPeriods)
            {
                //计算调整区间
                var modifiedPeriod =  deltaPeiod * 0.5;
                //防止进入提高效率
                if (period + modifiedPeriod < nextGuassPeriod)
                    continue;
                //计算频率
                var wn = SpectrumMethod.Getwn(period);
                //获得地震动反应谱对象
                var spectrumParam = this.waveSpectrum.GetSpectrumParam(period);
                //获得规范反应谱对应周期的对象
                var codeSpectrumParam = this.codeSpectrum.GetSpectrumParam(period);
                //不存在则继续
                if (codeSpectrumParam == null) continue;
                //计算相对误差判断是否满足
                var Eerror = Math.Abs((spectrumParam.Afa - codeSpectrumParam.Afa) / codeSpectrumParam.Afa);
                //计算误差函数
                var R_Error = codeSpectrumParam.Afa / spectrumParam.Afa;
                //如果误差满足则继续
                if (Eerror < ArtiWaveParam.ModifiedTol) continue;
                //获得响应
                var responseGlobal = this.waveSpectrum.GetResponseSet(spectrumParam.Period);
                //获得最大加速度响应值 及其对应的时刻
                var responseMax = responseGlobal.GetAfaMaxWithSignal(true);
                var responseMaxTime = responseGlobal.GetAfaMaxTime(true);
                //命中多少频率
                var testIndex = 0;
                var testList = new List<double>();
                //遍历
                for (; GuassIndex < this.StationaryGaussianList.Count; GuassIndex++)
                {
                    //不能为零
                    if (GuassIndex < 0) continue;
                    //获得对象
                    var stationaryGaussian = this.StationaryGaussianList[GuassIndex];
                    //如果周期比它大则停止
                    if (stationaryGaussian.period > period + modifiedPeriod)
                    {
                        nextGuassPeriod = stationaryGaussian.period;
                        GuassIndex--;
                        break;
                    }
                    else if (stationaryGaussian.period > period - modifiedPeriod)
                    {
                        //计算响应
                        var responseSub = stationaryGaussian.GetResponseSet(this.Ft, this.TimeInterval, responseMaxTime + 0.1, period);
                        //获得响应
                        var responseSubMax = responseSub.GetResponseParam(responseMaxTime).U * Math.Pow(wn, 2);
                        //设定贡献率
                        stationaryGaussian.ContributionFactor = Math.Abs(responseSubMax / responseMax);
                        //设定贡献率边界值
                        var contributionLimit = period <= ArtiWaveParam.Tg ? 
                            ArtiWaveParam.ShortPeriodContributionModified * ArtiWaveParam.ModifiedTol : 
                            ArtiWaveParam.ModifiedTol;
                        //满足贡献率
                        if (stationaryGaussian.ContributionFactor < contributionLimit) continue;
                        //获得符号
                        int sign = responseSubMax * responseMax < 0 ? -1 : 1;
                        //设定调整系数
                        stationaryGaussian.ModifyPowerSpectrum((float)Math.Pow(R_Error, sign));
                        testIndex++;
                    }
                }
            }
        }

        #region 赋值调整相位调整弃用方法

        /// <summary>
        /// 设定高斯平稳过程
        /// </summary>
        //private void SetStationaryGaussianBySamePeriod()
        //{
        //    //初始化
        //    this.StationaryGaussianList = new List<StationaryGaussianParam>();
        //    //构造随机数
        //    Random rd = new Random();
        //    //确定功率个数
        //    var num = (int)(ArtiWaveParam.PeriodLength / ArtiWaveParam.PeriodDelta);
        //    //遍历构造高斯平稳过程
        //    for (int index = 1; index <= num; index++)
        //    {
        //        //获得此可周期值
        //        var period = ArtiWaveParam.PeriodDelta * index;
        //        //求得功率增量
        //        var deltaW = 2 * Math.PI / period - 2 * Math.PI / (period + ArtiWaveParam.PeriodDelta);
        //        //获得基础相位角
        //        var testThetaK = 2 * Math.PI / 100 * rd.Next(100);
        //        //添加平稳高斯过程
        //        this.StationaryGaussianList.Add(new StationaryGaussianParam(
        //             2 * Math.PI / period, deltaW, testThetaK, this.Duration, this.codeSpectrum));
        //        //暂停
        //        Thread.Sleep(1);
        //    }
        //}

        /// <summary>
        /// 修改功率谱值
        /// </summary>
        //private void ModifyPowerSpectrum()
        //{
        //    //各控制周期的调整系数字典 R E
        //    var powerModifyDict = new SortedDictionary<double, Tuple<double, double>>();
        //    //遍历地震动反应谱
        //    for (int index = 0; index < this.waveSpectrum.Count; index++)
        //    {
        //        //获得地震动反应谱对象
        //        var spectrumParam = this.waveSpectrum[index];
        //        //获得规范反应谱对应周期的对象
        //        var codeSpectrumParam = this.codeSpectrum.GetSpectrumParam(spectrumParam.Period);
        //        //不存在则继续
        //        if (codeSpectrumParam == null) continue;
        //        //获得响应
        //        //var responseGlobal = this.waveSpectrum.GetResponseSet(spectrumParam.Period);
        //        //获得平稳高斯参数
        //        //var stationaryGaussian = this.GetStationaryGaussianParam(spectrumParam.Period);
        //        //不存在则继续
        //        //if (codeSpectrumParam == null || responseGlobal == null || stationaryGaussian == null) continue;
        //        //获得激励下响应
        //        //var responseSub = stationaryGaussian.GetResponseSet(this.Ft, this.TimeInterval, this.Duration, spectrumParam.Period);
        //        //获得符号（有问题！！！）
        //        //int sign = responseSub.GetAfaMaxWithSignal(true) * responseGlobal.GetAfaMaxWithSignal(true) < 0 ? -1 : 1;
        //        int sign = 1;
        //        //计算误差函数R
        //        var valueTuple = Tuple.Create<double, double>(
        //            Math.Pow((codeSpectrumParam.Afa) / spectrumParam.Afa, sign),
        //            (spectrumParam.Afa - codeSpectrumParam.Afa) / codeSpectrumParam.Afa);
        //        //计算相对误差
        //        powerModifyDict.Add(spectrumParam.Period, valueTuple);
        //    }
        //    //修正谱密度
        //    this.StationaryGaussianList.ForEach(gauss => gauss.ModifyPowerSpectrum(powerModifyDict, ArtiWaveParam.ModifiedTol));
        //}

        /// <summary>
        /// 仅调整某目标频率下的分量
        /// </summary>
        //private void ModifyPowerSpectrumTest3()
        //{
        //    //周期增量
        //    var deltaPeiod = (float)ArtiWaveParam.ModifiedPeriodDelta;
        //    //遍历地震动反应谱
        //    for (float period = deltaPeiod; period <= ArtiWaveParam.Tg; period += deltaPeiod)
        //    {
        //        //计算频率
        //        var wn = SpectrumMethod.Getwn(period);
        //        //获得地震动反应谱对象
        //        var spectrumParam = this.waveSpectrum.GetSpectrumParam(period);
        //        //获得规范反应谱对应周期的对象
        //        var codeSpectrumParam = this.codeSpectrum.GetSpectrumParam(period);
        //        //不存在则继续
        //        if (codeSpectrumParam == null) continue;
        //        //计算相对误差判断是否满足
        //        var Eerror = Math.Abs((spectrumParam.Afa - codeSpectrumParam.Afa) / codeSpectrumParam.Afa);
        //        //计算误差函数
        //        var R_Error = codeSpectrumParam.Afa / spectrumParam.Afa;
        //        //如果误差满足则继续
        //        if (Eerror < 4 * ArtiWaveParam.ModifiedTol) continue;
        //        //获得响应
        //        var responseGlobal = this.waveSpectrum.GetResponseSet(spectrumParam.Period);
        //        //获得最大加速度响应值 及其对应的时刻
        //        var responseMax = responseGlobal.GetAfaMaxWithSignal(true);
        //        var responseMaxTime = responseGlobal.GetAfaMaxTime(true);
        //        //满足贡献率的列表（测试）
        //        var posList = new List<StationaryGaussianParam>();
        //        var negList = new List<StationaryGaussianParam>();
        //        //不满足则求解各功率的时程响应 判断是否贡献较大
        //        this.StationaryGaussianList.ForEach(stationaryGaussian =>
        //        {
        //            //忽略小分量
        //            if (stationaryGaussian.period > period - deltaPeiod * 0.1 && stationaryGaussian.period < period + deltaPeiod * 0.1)
        //            {
        //                //计算响应
        //                var responseSub = stationaryGaussian.GetResponseSet(this.Ft, this.TimeInterval, responseMaxTime + 1, period);
        //                //获得响应
        //                var responseSubMax = responseSub.GetResponseParam(responseMaxTime).U * Math.Pow(wn, 2);
        //                //如果响应过小调整相位
        //                if (Math.Abs(responseSubMax / responseMax) < ArtiWaveParam.ContributionTol)
        //                {
        //                    stationaryGaussian.thetak += Math.PI;
        //                }
        //            }
        //        });
        //        //不满足则求解各功率的时程响应 判断是否贡献较大
        //        //this.StationaryGaussianList.ForEach(stationaryGaussian =>
        //        //{
        //        //    stationaryGaussian.Add(this.Ft, responseMax, responseMaxTime, this.TimeInterval, period, Eerror, R_Error);
        //        //});
        //    }
        //    //this.StationaryGaussianList.ForEach(stationaryGaussian => stationaryGaussian.ModifyPowerSpectrum());
        //}

        //private void ModifyPowerSpectrumTest4(float period)
        //{

        //    //计算频率
        //    var wn = SpectrumMethod.Getwn(period);
        //    //获得地震动反应谱对象
        //    var spectrumParam = this.waveSpectrum.GetSpectrumParam(period);
        //    //获得规范反应谱对应周期的对象
        //    var codeSpectrumParam = this.codeSpectrum.GetSpectrumParam(period);
        //    //不存在则继续
        //    if (codeSpectrumParam == null) return;
        //    //计算相对误差判断是否满足
        //    var Eerror = Math.Abs((spectrumParam.Afa - codeSpectrumParam.Afa) / codeSpectrumParam.Afa);
        //    //计算误差函数
        //    var R_Error = codeSpectrumParam.Afa / spectrumParam.Afa;
        //    //如果误差满足则继续
        //    if (Eerror < ArtiWaveParam.ModifiedTol) return;
        //    //获得响应
        //    var responseGlobal = this.waveSpectrum.GetResponseSet(spectrumParam.Period);
        //    //获得最大加速度响应值 及其对应的时刻
        //    var responseMax = responseGlobal.GetAfaMaxWithSignal(true);
        //    var responseMaxTime = responseGlobal.GetAfaMaxTime(true);
        //    //满足贡献率的列表（测试）
        //    var posList = new List<StationaryGaussianParam>();
        //    var negList = new List<StationaryGaussianParam>();
        //    //不满足则求解各功率的时程响应 判断是否贡献较大
        //    this.StationaryGaussianList.ForEach(stationaryGaussian =>
        //    {
        //        stationaryGaussian.Add(this.Ft, responseMax, responseMaxTime, this.TimeInterval, period, Eerror, R_Error);
        //    });
        //    this.StationaryGaussianList.ForEach(stationaryGaussian => stationaryGaussian.ModifyPowerSpectrum());
        //}

        /// <summary>
        /// 获得平稳高斯参数
        /// </summary>
        /// <param name="period"></param>
        //private StationaryGaussianParam GetStationaryGaussianParam(float period)
        //{
        //    //功率
        //    //var w = Math.PI * 2 / period;
        //    //根据周期寻找最近的谱值
        //    var responses = this.StationaryGaussianList.Find("period", period, ArtiWaveParam.PeriodDelta / 10);
        //    //是否存在多个值
        //    if (responses.Count == 0)
        //    {
        //        //MessageBoxExtern.Error(string.Format("目标谱周期{0}对应谱值不存在", period));
        //        return null;
        //    }
        //    else if (responses.Count > 1)
        //        MessageBoxExtern.Error(string.Format("目标谱周期{0}存在相似谱值", period));
        //    return responses.First();
        //}

        #endregion

        /// <summary>
        /// 绘图
        /// </summary>
        /// <param name="waveChart"></param>
        /// <param name="spectrmChart"></param>
        public void ToChart(Chart waveChart, Chart spectrmChart)
        {
            this.accels.ToChart(waveChart);
            this.waveSpectrum.ToChart(spectrmChart, false);
            this.codeSpectrum.ToChart(spectrmChart, true);
        }

        /// <summary>
        /// 迭代
        /// </summary>
        private void IterationRun()
        {
            //最大迭代次数
            var IterationMax = (ArtiWaveParam.ThetaModifiedIntervalTime * ArtiWaveParam.IterationGroupNum);
            //首次迭代不考虑相位调整
            if (this.IterationCount != 0)
            {
                ////迭代次数
                //this.IterationCount++;
                ////调整相位
                //this.ModifyPowerSpectrumTest3();
                ////this.ModifyPowerSpectrumTest4(0.1F);
                ////初始化非平稳随机规程
                //this.SetAccels();
                ////完成反应谱求解
                //this.SetSpectrum();
            }
            //三次幅值谱调整
            for (int i = 1; i <= ArtiWaveParam.ThetaModifiedIntervalTime; i++)
            {
                //迭代次数
                this.IterationCount++;
                //修改幅值谱
                this.ModifyPowerSpectrumA();
                //初始化非平稳随机规程
                this.SetAccels();
                //完成反应谱求解
                this.SetSpectrum(IterationMax == this.IterationCount);
            }
        }

        /// <summary>
        /// 迭代测试
        /// </summary>
        /// <param name="cout"></param>
        public void IterationGroupRun(int progressIndex, Action<int, float> setProgress)
        {
            for (int index = 0; index < ArtiWaveParam.IterationGroupNum; index++)
            {
                IterationRun();
                setProgress(progressIndex,
                    (float)this.IterationCount / (ArtiWaveParam.ThetaModifiedIntervalTime * ArtiWaveParam.IterationGroupNum));
            }
            //释放
            this.waveSpectrum.TrimExcess();
            this.StationaryGaussianList.Clear();
            this.StationaryGaussianList.TrimExcess();
            GC.Collect();
        }

        /// <summary>
        /// 构造函数
        /// </summary>
        public ArtiWaveObj(SpectSet spectrum, string name, int tmin, int tmax, int duration, double timeInterval)
        {
            //获得人工波名称
            this.Name = name;
            //获得持时与时间间隔
            this.Duration = duration; this.TimeInterval = timeInterval; this.TimeIndexMax = (int)(this.Duration / this.TimeInterval);
            //设定规范反应谱
            this.codeSpectrum = spectrum;
            //设定包络函数
            this.Ft = new EnvelopeFunction(tmin, tmax, ArtiWaveParam.c);
            //设定高斯平稳过程
            this.SetStationaryGaussianBySameW();
            //初始化非平稳随机规程
            this.SetAccels();
            //完成第一次反应谱求解
            this.SetSpectrum(false);
        }

        /// <summary>
        /// 空的构造函数
        /// </summary>
        public ArtiWaveObj()
        {

        }
    }
}
